package ml;



public class MicrolensingSimulation {
    static public final double sigmaStars = 0.064;
    static public final double sigmaMatter = 0.576;
    static public final double shear = 0.52;

    //    static public final double sigmaStars = 0.0;
    //    static public final double sigmaMatter = 0.0;
    //    static public final double shear = 0.0;

    static public final double sigmaTotal = sigmaStars+sigmaMatter;

    //Map size
    static private double sizeMap = 16.0;
    //    static private final double sizeMap = 8.0;
    static private final int nPixMap = 1024;
    static private final double pixSize = sizeMap/nPixMap;

    static private final double   margin = 1.25;

    static public  final int  nRays = 6; //So there is aprox. 1000 rays/pixel
    static private int numberOfProcessors = Runtime.getRuntime().availableProcessors();
    //  static private int numberOfProcessors = 1;

    private static void initArray(double[][] array){
	int sizex = array[0].length;
	int sizey = array.length;

	System.out.println("initializing array sizex: " + sizex + ",sizey: " + sizey);
	for (int j=0;j<sizey;j++){
	    for (int i=0;i<sizex;i++){
		array[j][i] = 0;
	    }
	}

    }

    /**
     * @param args
     */
    public static void main(String[] args) {
	double[][] magnificationMap = new double[nPixMap][nPixMap];
	initArray(magnificationMap);
	// As we need a square map we are going to shoot rays in
	//The following region
	double x1 =sizeMap/Math.abs(1.-sigmaTotal-shear);
	double x2 =sizeMap/Math.abs(1.-sigmaTotal+shear);
	double sizeX =margin* Math.max(x1, x2);
	double shootRegionAspectRatio = Math.abs(1.-sigmaTotal+shear)/Math.abs(1.-sigmaTotal-shear);
	double sizeY = sizeX * shootRegionAspectRatio;


	int shootRegionX=(int)Math.round(margin * nPixMap);
	int shootRegionY=(int)Math.ceil(shootRegionX* shootRegionAspectRatio);

	double[][] shootingMap = new double[shootRegionY][shootRegionX];
	initArray(shootingMap);
	System.out.println("Number of cores: " + numberOfProcessors);
	System.out.println("Size x: " + sizeX);
	System.out.println("Size y: " + sizeY);
	double stepX = pixSize/nRays;
	double stepY = pixSize*shootRegionAspectRatio/nRays;
	System.out.println("Step X: " + stepX);
	System.out.println("Step Y: " + stepY);

	//Number of stars
	int nStars = (int)Math.round(sigmaStars*(sizeX * sizeY));
	System.out.println("Numero de estrellas: " + nStars);
	double xStar[]={-sizeMap/4,0.0, +sizeMap/4,0.0 };
	double yStar[]={-sizeMap*Math.tan(Math.PI/6)/4,sizeMap/(Math.cos(Math.PI/6)*4),-sizeMap*Math.tan(Math.PI/6)/4,0.0};
	StarBag starBag = new StarBag(nStars, sizeX, sizeY);
	//	StarBag starBag = new StarBag(xStar,yStar);

	long initialTime = System.currentTimeMillis();
	//		double progress, elapsedTime, remaining;

	int ny = (int)Math.round(sizeY/stepY);
	sizeY = ny *stepY;
	int nyChunk = ny/numberOfProcessors;
	Deflector[] deflectors = new  Deflector[numberOfProcessors];

	for (int i = 0;i<numberOfProcessors;i++){
	    double startY=-sizeY/2 + i*nyChunk*stepY;
	    double endY = startY + nyChunk*stepY;
	    deflectors[i] = new Deflector(magnificationMap, shootingMap, starBag, 
		    -sizeX/2, sizeX/2, stepX,
		    startY, endY, stepY,sizeX, sizeY, nPixMap, pixSize, sigmaMatter, shear);
	}

	//Start all the deflectors
	for (int i = 0;i<numberOfProcessors;i++){
	    deflectors[i].start();
	}	

	boolean working = true;
	//Wait until all have been finished;
	while (working){
	    try {
		Thread.sleep(10000);
	    } catch (InterruptedException e) {
		// TODO Auto-generated catch block
		e.printStackTrace();
	    }
	    working = false;
	    for (int i = 0;i<numberOfProcessors;i++){
		working |= deflectors[i].isAlive();
		System.out.println(deflectors[i].toString());
	    }	
	}

	//Compute flux parameters
	long totalFlux =0;
	long lostFlux =0;

	for (int t = 0; t<numberOfProcessors;t++){
	    totalFlux += deflectors[t].totalFlux();
	    lostFlux += deflectors[t].lostFlux();
	}

	System.out.println("Total flux: " + totalFlux);
	System.out.println("Lost flux: " + lostFlux);
	System.out.println("Lost fraction: " + ((double)lostFlux)/totalFlux);


	//Calculate magnification
	double meanMagnification =0;
	for (int j=0;j<nPixMap;j++){
	    for (int i=0;i<nPixMap;i++){
		magnificationMap[j][i] /=(nRays*nRays);
		meanMagnification += magnificationMap[j][i];
	    }

	}
	meanMagnification/=(nPixMap*nPixMap);
	System.out.println("Mean magnification: " + meanMagnification);
	//Saving FITS
	String[] keys ={"SIG_STAR", "SIG_DIST", "SHEAR", "SIZE", "MARGIN", "NRAYS", "NSTARS"};
	String[] values ={Double.toString(sigmaStars), 
		Double.toString(sigmaMatter), 
		Double.toString(shear),
		Double.toString(sizeMap), 
		Double.toString(margin), 
		Integer.toString(nRays*nRays),
		Integer.toString(starBag.size())};
	String[] comments ={"Star density", 
		"Extended mass density", 
		"Shear parameter",
		"Map size in Einstein radii", 
		"Shoot region margin factor", 
		"Number of rays per pixel",
	"Number of stars in shooting region"};
	Utils.saveFits(magnificationMap,shootingMap,starBag, "magPattern.fits", keys, values, comments);
	//	Utils.saveFits(, "shootingPattern.fits", keys, values, comments);


    }

}
